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Abstract 

The aim of this paper is to give an explicit formula of the invariant distribution of a 
quasi-birth-and-death process in terms of the block entries of the transition probability 
matrix using a matrix-valued orthogonal polynomials approach. We will show that the 
invariant distribution can be computed using the squared norms of the corresponding 
matrix- valued orthogonal polynomials, no matter if they are or not diagonal matrices. We 
will give an example where the squared norms are not diagonal matrices, but nevertheless 
we can compute its invariant distribution. 

1 Introduction 

The connection between random walks/birth-and-death processes and orthogonal polynomi- 
als is very well known. In both cases the state space is the set of nonnegative integers, i.e. 
S = {0, 1, . . .}, while the parameter set is T = {0, 1, . . .} for random walks (i.e. discrete time) 
and T = [0,oo) for birth-and-death processes (i.e. continuous time). In a series of papers, 
Karlin and McGregor (see [101 [11]) found an appropriate tool to study these processes 
connecting the tridiagonal one-step transition probability matrix P (or the infinitesimal gen- 
erator A in the continuous time case) with a measure supported in the real line. In particular, 
they obtained an integral representation for the n-step transition probability matrix P n (or 
the transition probability matrix P{t) in the continuous case) in terms of the spectral mea- 
sure and the corresponding orthogonal polynomials, as well as the explicit expression of the 
invariant distribution, i.e. a row vector tv with nonnegative components such that nP = ir. 

Quasi-birth-and-death processes are a natural extension where now the state space is two 
dimensional, i.e. S = % > 0, 1 < j < N}. The first component of the pair is usually 
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called the level and the second one the phase. Now the one-step transition probability matrix 
P is block tridiagonal. The link with matrix-valued orthogonal polynomials was initially 
raised independently by [2] and [3] in discrete time (for continuous time see [3]), although 
the first example may be traced back to the last section of where the authors deal with 
the case of a random walk with state space the set of all integers, replacing the spectral 
measure by a 2 x 2 non-negative matrix. In [2] and [1] one can find an integral representation 
for the n-step transition probability matrix P n , along with other probabilistic useful results. 
For a much more detailed presentation of this field, as well as its connections with queueing 
problems in network theory and the general field of communication systems the reader should 
consult P31 HUE]. 

Nevertheless computing the invariant distribution of a quasi-birth-and-death process, i.e. 
a row vector 7r with nonnegative components such that tvP = 7T, is a much harder problem, 
compared with the 1-dimensional situation. It is possible to derive nicer looking expressions 
for the invariant distribution for some special cases, like level-independent quasi-birth-and- 
death processes (see [12] ) . Also in some papers (see [5j [6J El E]) the invariant distribution 
was given in terms of the entries of the diagonal norms of the corresponding matrix-valued 
orthogonal polynomials, giving a natural candidate. Here we will show that it is possible 
to derive an explicit expression of the invariant distribution in terms of the block entries 
of the transition probability matrix or equivalently, in terms of the squared norms of the 
corresponding matrix-valued orthogonal polynomials, no matter if they are or not diagonal. 
In fact, we will give an example of a quasi-birth-and-death process where the squared norms 
of the corresponding matrix- valued orthogonal polynomials are not diagonal, but nevertheless 
we can compute its invariant distribution. 

In Section [2] we will define the matrix version of the so-called potential coefficients and 
connect them with the squared norms of matrix-valued orthogonal polynomials. In Section 
[3] we will give the explicit expression of an invariant distribution in terms of these matrix- 
valued potential coefficients and finally, in Section [U we will give the example mentioned in 
the paragraph above. 

2 Matrix-valued potential coefficients 

The content in this section is already known in the literature, but we will point out some 
properties that these matrix-valued potential coefficients have in order to prove the main 
result in the next section. For simplicity, we will focus on the case of discrete time quasi- 
birth-and-death processes. The continuous time case only suffers cosmetic changes. 

Consider a nonhomogeneous discrete time quasi-birth-and-death process, i.e. a discrete 
time Markov process on the countable two dimensional state space S = : i > 0,1 < 

j < N} with block tridiagonal transition probability matrix 




P = 



(B A 
Ci B 1 Ai 

C2 B>2 A2 



V 
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The first component of the pair is usually called the level and the second one the phase. P 
is stochastic, i.e. all entries are nonnegative and all rows sum up to one, that is (Bq + A^je^ = 
ejy and (C n + B n + A n )e^ = ejy, n > 1, where e^r is the column vector of l's of dimension 
N. We will assume that our process is irreducible and that the coefficients A n and C n are 
nonsingular matrices. Clearly in the case when the number of phases N is one we are back 
to the case of an ordinary birth-and-death process. 

In order to link these processes with matrix-valued orthogonal polynomials the matrix P 
needs to be transformed into a symmetric matrix. This is possible if there exists a nonsingular 
block diagonal matrix 

Rq 

r= I Ri 



such that RPR 1 is symmetric (see for instance Theorem 2.1 of [2]). The matrices R n are 
subject to the following restrictions 

(2.2) RnBnR^ 1 = {R n B n R n l ) J ', and R n A n R n ^ 1 = {R n+ iC n+ iR n 1 ) T , n > 0, 

where M T denotes the transpose of the matrix M. 
Let us call 

(2.3) II n = R^Rn, n > 0, 

which are obviously symmetric. From (12. 2j) we have that n n , n > 0, satisfy 

(2.4) n n B n = B?n n , n > 0, 
and 

(2.5) U n A n = Cj +1 II n+lj n>0. 

Condition (|2.5p gives an explicit formula for Tl n , n > 1, given A n , C n and LTo: 

(2.6) n n = (C 1 T C 2 T ---Cj)- 1 n (A A 1 ---A n _ 1 ), n>l, 

This formula is similar to the scalar potential coefficients for birth-and-death processes (see 
[10|). so we will call them matrix-valued potential coefficients. 

Under these assumptions, there always exists a weight matrix W such that the matrix- 
valued polynomials defined by the three-term recurrence relation 

(2.7) xQ n (x) = A n Q n+1 (x) + B n Q n (x) + C n Q n -x(x), n > 0, 
where Q-i(x) = and Qo(x) = I, are orthogonal, i.e. 

Q n (x)W(x)Q^(x)dx = \\Qn\\w&nm- 
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In particular, as it is remarked in page 121 of [2], IIo can be given in terms of the 0-th moment 
So = fW(x)dx. For an orthonormal family (Q n ) n we have that 

1= j Qo{x)W{x)Ql{x)dx = R S Rl 

Therefore 

n = 5 - 1 = (||Q ||^)- 1 . 

That means that for every weight matrix W associated with the Jacobi matrix (12. ip . the 
matrix-valued potential coefficients are defined recursively using (|2.6p in a unique way. 

Remark 2.1. Observe that the sequence (R n )n is not unique in the sense that we can always 
consider any other sequence (U n R n ) n , where (U n ) n is any sequence of orthogonal matrices 
and we will get another family of orthonormal polynomials. Nevertheless, the matrix-valued 
potential coefficients II n are always unique (see (12.3P ). 

Finally we will show that the matrix-valued potential coefficients can be given in terms of 
the squared norms of the matrix-valued orthogonal polynomials. From the three-term recur- 
rence relation (|2.7p we can multiply on the right by WQ n , WQ n+ i and WQ n -i, respectively, 
and integrate in such a way that the sequence of squared norms satisfies the following 

relations 

-SnllQnllvy = IIQn.||u/^n) II Qn \\ = || Qn+1 \\ W i n > 0. 

But these are just the relations (12. 4p and (12. 5p . Since IIo = (IIQoIIvf) -1 ; we have that 
(2.8) n n = {\\Q n \\w)'\ n>0. 

In particular we have that the matrix-valued potential coefficients II n are positive semi- 
definite matrices. 



3 The invariant distribution 



With all the information from the last section we are able to give a very natural candidate 
for the invariant distribution of the process P using the matrix- valued potential coefficients 

n n . 

Theorem 3.1. Let P be the transition probability matrix given by (12. ip . Define the sequence 
of matrices H n , n > 1, as in (|2.6p with Ho = (J* W(x)dx)~ 1 or, equivalently, as in (|2.8p . 
Consider the following row vector 

(3.1) TV = ((U e N ) T ; (U ieN f; (U 2 e N f; • • • ), 

where is the column vector of 1's of dimension N. Then n is an invariant distribution 
for the process P, i.e. all components of it are nonnegative and 

(3.2) 7rP = 7T. 
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Proof. All components of 7r are nonnegative since Il n are positive semi-definite matrices (see 

To prove (|3.2p . we have to check that 

(U e N ) T B + (IliejvfCi = (n ejv) T 

and 

{Ii n -ie N ) T A n ^i + (Il n e N ) T B n + (Il n+1 e N ) T C n+1 = (U n e N ) T , n > 1. 
The first equality holds using properties (|2.4p . ()2,5p and the fact that P is stochastic, since 

e^(rioA) + n x Ci) = e^(5 T n + ^n ) = [(B + A )e N ] T U = (U e N ) T , 
Similarly, for n > 1 

e^(n„_ 1 A n _ 1 + n„5„ + n„ +1 c„, + i) = [(A n + B n + C n )e N ] T U n = {U n e N ) T . 

□ 

Remark 3.1. The same result holds in the continuous time case. In this case, the transition 
probability matrix P is replaced by an infinitesimal generator A, block tridiagonal as in (|2.ip . 
but with (Bq + Ao)ejv = and (C n + B n + A n )e^ = 0, n > 1 ('see J3j/ /or more details). 
Therefore the row vector n defined by (|3.ip is an invariant distribution for the process, i.e. 
all components are nonnegative and 

tvA = 0. 

Let us make a few comments about the unicity of this invariant distribution. If the process 
is recurrent then there exists a unique invariant distribution given by ()3.ip (see Theorem 5.4 
of |16j). But if the process is transient there are no results about unicity. In fact, the case 
of a random walk on the integers treated in [TT] gives (for p ^ q) an example of a transient 
process where the invariant distribution is not unique. 

In Section [2] we saw that the matrix- valued potential coefficients U n are defined in terms 
of the corresponding weight matrix W and the block entries of P. These quantities are 
constructed in a unique way, so the question of the unicity of the invariant distribution 
is the same as the question of the unicity of the weight matrix W corresponding to the 
equivalent class of block tridiagonal Jacobi matrices P by unitary transformations. So if we 
can guarantee that the corresponding weight matrix is unique, then the invariant distribution 
will be also unique. Nevertheless there are not general results to determine if the weight 
matrix associated with a block tridiagonal Jacobi matrix is unique. 

Observe that Theorem 13.11 gives a way to compute an invariant distribution even if the 
squared norms of the corresponding matrix-valued orthogonal polynomials are not diagonal. 
Until now, all the examples introduced in the literature had the special lucky property that 
these norms were diagonal matrices. In the following section we will give an example of a 
quasi-birth-and-death process where the squared norms of the corresponding matrix-valued 
orthogonal polynomials are not diagonal, but nevertheless, using (|3.1[) . we can compute its 
invariant distribution. This example, as far as the author knows, may be the first quasi- 
birth-and-death process with this property in the literature. 
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4 The example 



Consider the following 2x2 weight matrix 

w( \ a/i ,R(kx 2 + P-k + l {p-k + l){l-x)\ rn 

where a,/3 > — 1 and < k < /3 + 1. This weight matrix was introduced in [Tj in a different 
context, that of searching second-order differential operators having several families of matrix- 
valued orthogonal polynomials as eigenfunctions. Let us call (P n )n the unique monic family 
of matrix-valued orthogonal polynomials associated with W. 
Let Ao be the nonsingular matrix 

/-, oc+P-k+3 \ 
A _ I a+2/3-2fc+4 ] 



1 



«+2/3-2fc+4 I ! 
P-k+l J 



and let us consider the equivalent weight matrix 

W = A WA%. 



Now we consider a special family of matrix-valued orthogonal polynomials with respect to 
W of the following form 

Q n (x) = A n P n (x)AQ 1 , 

where 

,2 



A, 



a n (n 2 + n(a + g + 3) + k(a + 2/3 - 2k + 4)) 

k(a + f3-k + n + 3) "' 



with 



, b n {n 2 + n(a + g + 3) + fc(a + 2/3 - 2fc + 4 )) 

(fe + n)(/3-fc + l) 

(a+/3+n+3) n (a+/3-fc+n+3)(n(a+/3+n+2)+fc(q+l)) 

° n ~~ (/3+2)„(n(a+/3) 2 +n(2n+5)(a+/3)+n(ri+2)(n+3)+fc(2a^+2/3-fc(ri+2)+a 2 +5a-2fca-n 2 ^ ' 

b (a + g + n + 3) n (fc + n)(n(a + g + n + 2) + fc(a + 1)) 

n ~ [fi + 2) n ((n 2 + 2nk)(a + (3) + n(n + 5)k + k 2 (a - n + 1) + n 2 (n + 3)) ' 

and (z) n will denote the Pochhammer symbol defined by (z) n = z(z + 1) ■ ■ ■ (z + n — 1) for 
n > 0, (z)o = 1. Observe that Qo = I- The choice of the leading coefficient A„Aq 1 of the 
family of polynomials (Q n )n is motivated by the remarkable fact that 

(4.1) Qn(l)e 2 = e 2 , 

where e 2 is the 2-dimensional column vector with all entries equal to 1. In other words, the 
sum of the elements in each row of Q n (l) gives the value 1. 
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The family of matrix-valued orthogonal polynomials introduced above satisfies a three- 
term recursion relation 



(4.2) xQ n (x) = A n Q n+ i(x) + B n Q n (x) + C n Q n -i{x), n > 0, 

where Q-i(x) = and Qo(x) = I and A n ,B n ,n > 0, C n ,n > 1 are 2x2 full matrices with 
nonnegative entries. The corresponding Jacobi matrix 



(4.3) 



P 



/B A \ 

d B x A x 

C 2 B2 A2 

v •• '• '•/ 



is stochastic. This is a consequence of choosing the family {Q n ) n with the property (|4.ip . 
Indeed, applying both sides of the identity (14, 2p to the vector e 2 , setting x = 1 and using 
(|4.ip we obtain that (Bq + ^4o) e 2 = e 2 and (C n + B n + ^4 n )e2 = e2, n > 1, i.e. the sum of 
the entries in each row of P equals one. Therefore, it gives a quasi-birth-and-death process 
with two phases (N = 2) and depending on three parameters, a, (3 and k. 

The state space and the corresponding one-step transitions look as follows (see [H [8] 
for an explanation of the labeling and ordering of the states) 




Finally the squared norms ||Qn|||^ are 2x2 full matrices (not diagonal). Therefore, using 
Theorem 13.11 we get that an invariant distribution for this process is given by 

(4.4) TV = ((n e 2 ) T ; (n ie2 ) T ; (n 2 e 2 ) T ; • • • ), 

where 

(4-5) n n = (\\QnW~r 1 , n>0. 

We can use Theorem 8.1 of [8J to study the recurrence of the process, since our weight matrix 
has similar properties to that introduced in [8]. The Markov process that results from P is 
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never positive recurrent. If —1 < f3 < then the process is null recurrent. If /? > then 
the process is transient. For — 1 < f3 < the invariant distribution will be unique, as we 
mentioned at the end of Section [3j 

The explicit expressions of A n , B n , C n and ||Qn||~ are too long to be displayed here. 

Instead, we will give all these data for the special case of a = (3 = and k = \. Similar 
expressions can be derived for other values of the parameters a, f3 and k. The coefficients of 
the three-term recurrence relation (14.2p are, for n > 

/ (n + 3) 2 (2n 2 + 4n + l)(4n 2 + 18n + 17)(16n 4 + 128n 3 + 348n 2 + 368n + 117) 

A 2(n + 2)(2n + 3)(2n 2 + 8n + 7)(4n 2 + lOre + 3)(16n 4 + 160n 3 + 572n 2 + 860n + 463) 

2(n + 3) 2 (2n 2 + 4n + l) 2 (4n 2 + 18n + 17) 

V (2n + 3)(2n 2 + 8n + 7)(4n 3 + 14n 2 + 9n + l)(16n 4 + 160n 3 + 572n 2 + 860n + 463) 

2(n + 3) (2n 2 + 4n + 1) (2n 2 + 12n + 17) (4n 3 + 26n 2 + 49n + 28) \ 

(n + 2)(2n + 3)(2n 2 + 8n + 7)(4n 2 + lOre + 3)(16n 4 + 160n 3 + 572n 2 + 860n + 463) 
(n + 3)(2n 2 + An + l)(4n 3 + 26n 2 + 49n + 28)(16n 4 + 128n 3 + 348n 2 + 368n + 125) 

2(2n + 3)(2n 2 + 8n + 7)(4n 3 + 14n 2 + 9n + l)(16n 4 + 160n 3 + 572ra 2 + 860n + 463) / 

/ 32n 8 + 384n 7 + 1928n 6 + 5256n 5 + 8450n 4 + 8148n 3 + 4577n 2 + 1365n + 163 

R = (2n 2 + 4n + l)(2n 2 + 8n + 7)(16n 4 + 96n 3 + 188n 2 + 132n + 31) 

(n + 2)(4n 2 + lOn + 3)(32n 6 + 256n 5 + 760n 4 + 1040n 3 + 674n 2 + 186n + 13) 

V 2(2n 2 + 4n + l)(2n 2 + 8n + 7)(4n 3 + 14n 2 + 9n + l)(16n 4 + 96n 3 + 188n 2 + 132n + 31) 

(4n 3 + 14n 2 + 9n + l)(32n 6 + 320n 5 + 1240n 4 + 2320n 3 + 2114n 2 + 834n + 121) 

2(n + 2)(2n 2 + 4n + l)(2n 2 + 8n + 7)(4n 2 + lOn + 3)(16n 4 + 96n 3 + 188n 2 + 132n + 31) 
(2n 2 + 6n + 3)(4n 3 + 18n 2 + 21n + 6)(4n 3 + 18n 2 + 21n + 3) 

(2n 2 + 4n + l)(2n 2 + 8n + 7)(16n 4 + 96n 3 + 188n 2 + 132n + 31) 

and for n > 1 

/ n(n + l)(4ra 2 + 2n - 3)(32n 6 + 192n 5 + 328n 4 + 32n 3 - 226n 2 - 4n + 33) 

r = 2(n + 2)(2n + 3)(2n 2 - l)(4n 2 + lOn + 3)(16n 4 + 32n 3 - 4n 2 - 20n + 7) 
n(n + l)(4n 2 + 2n - 3)(8n 4 + 32n 3 + 36n 2 + 8n - 5) 

V (2n + 3)(2n 2 - l)(4n 3 + 14n 2 + 9n + l)(16n 4 + 32n 3 - 4n 2 - 20n + 7) 

n(4n 3 + 2n 2 - 7n + 2) (8n 4 + 32n 3 + 36n 2 + 8n - 5) \ 

(n + 2)(2n + 3)(2n 2 - l)(4n 2 + lOn + 3)(16n 4 + 32n 3 - 4n 2 - 20n + 7) 
n(4n 3 + 2n 2 - 7n + 2)(32n fe + 192n 5 + 392n 4 + 288n 3 - 34n 2 - 132n - 43) 

2(2n + 3)(2n 2 - l)(4n 3 + 14n 2 + 9n + l)(16n 4 + 32n 3 - 4n 2 - 20n + 7) / 

We see that all entries of the transition probability matrix fj4.3f) are nonnegative rational 
functions in n, and it is very remarkable that the sum of the entries in each row of P equals 
one. 
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The corresponding squared norms are the 2x2 non-diagonal matrices 



WQr 



l~ 



/ (2n 2 + 4n + l)(16n 6 + 160n 5 + 628n 4 + 1212n 3 + 1173n 2 + 514n + 79) 

(n + l)(2n + 3)(4n 2 + lOn + 3) 2 (n + 2) 3 
(2n + l)(2n + 5)(2n 2 + 4n + 1) 

\ ~ 2(n + l)(n + 2) 2 (4n 2 + lOn + 3)(4n 3 + 14n 2 + 9n + 1) 

(2n + l)(2n + 5)(2n 2 + 4n + 1) 

~ 2(n + l)(n + 2) 2 (4n 2 + lOn + 3)(4n 3 + 14n 2 + 9n + 1) 
(2n 2 + 4n + l)(16n 6 + 128n 5 + 388n 4 + 564n 3 + 417n 2 + 152n + 22) 

(n + l)(n + 2)(2n + 3)(4n 3 + 14n 2 + 9n + l) 2 



/ 

The unique (the process is null recurrent in this case) invariant distribution is given by 
(jOl) where II n e2 can be calculated using (14,51) 



n„e 2 



/ 2(ra + l) 2 (n + 2) 2 (2n + 3)(4n 2 + 10n + 3)(4n 2 + 14n + 9) \ 

(2n 2 + 4n + l)(2n 2 + 8n + 7)(16n 4 + 96n 3 + 188n 2 + 132n + 31) 
2(n + l)(n + 2)(2n + 3)(4n 3 + 14n 2 + 9n + l)(4n 3 + 22n 2 + 33n + 8) 

V (2n 2 + 4n + l)(2n 2 + 8n + 7)(16n 4 + 96n 3 + 188n 2 + 132n + 31) / 



n > 0. 



If we denote by 7r™ and 77^ the two components of II n e2, we can study the behavior of 
the invariant distribution in Figure 4.1, a luxury we can afford since we have an analytic 
expression. 




1 2 3 4 5 6 
Figure 1: a = 0, /3 = 0, k = 0.5 
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